#### �p�b�P�[�W�̓ǂݍ��� ####
require(weights)
require(DirichletReg)
require(compositions)
require(MASS)
require(stringi)
require(ltm)
require(maxLik)
require(pbivnorm)
require(survey)

#### �f�[�^�̓ǂݍ��� ####
# ��N�w����
young.data.raw <- read.csv("young_data.csv")
# �S�N���
adult.data.raw <- read.csv("adult_data.csv")

## �s���Ӊ񓚎҂̏��O
young.data <- subset(young.data.raw, ! ((young.data.raw$age < 18 & young.data.raw$education > 3) | 
                                          (young.data.raw$age < 21 & young.data.raw$education == 7)))
adult.data <- subset(adult.data.raw, check.1 > 4 & check.2 < 4)

# �s���Ӊ񓚎҂̐l���Ɗ���
nrow(young.data.raw) - nrow(young.data)
round((nrow(young.data.raw) - nrow(young.data)) / nrow(young.data.raw), 3) * 100

sum(adult.data.raw$check.1 < 5)
round(sum(adult.data.raw$check.1 < 5) / nrow(adult.data.raw), 3) * 100
sum(adult.data.raw$check.2 > 3)
round(sum(adult.data.raw$check.2 > 3) / nrow(adult.data.raw), 3) * 100
(nrow(adult.data.raw) - nrow(adult.data))
round((nrow(adult.data.raw) - nrow(adult.data)) / nrow(adult.data.raw), 3) * 100

## ��v�ȕ��͂Ŏg��15���ڂ̑��_�ԓx�f�[�^�̒��o
# .raw�́u�킩��Ȃ��v�������l�������Ȃ��f�[�^
young.issue.data.raw <- young.issue.data <- young.data[, 1:15]
young.issue.data[young.issue.data == 6] <- NA
adult.issue.data.raw <- adult.issue.data <- adult.data[, 1:15]
adult.issue.data[adult.issue.data == 6] <- NA

issue.labels.en <- c("article.9", "defense", "CDF", "North.Korea", "apology", "nationality", 
                     "moral", "security", "death.penalty", "gender", "surnames", "family", "gay", 
                     "nuclear", "equality")
colnames(young.issue.data.raw) <- colnames(young.issue.data) <- 
  colnames(adult.issue.data.raw) <- colnames(adult.issue.data) <- issue.labels.en

## �����ƒ����f�[�^
UTAES.2017 <- read.csv("2017UTASP20180628.csv")

# �����l����
UTAES.2017[UTAES.2017 == 66 | UTAES.2017 == 99] <- NA

# ���I�҂݂̂𒊏o
elite.data <- subset(UTAES.2017, RESULT > 1)

# ���I�ғ��ł̉�����i�񓚎҂�RESPONSE = 1�j
round(prop.table(table(elite.data$RESPONSE)), 3)

#### �L�q���v ####
## ���_�ԓx�̕��z
young.issue.distribution <- t(apply(young.issue.data, 2, wpct, weight = young.data$weight))
adult.issue.distribution <- t(apply(adult.issue.data, 2, wpct, weight = adult.data$weight))

## �}1
issue.labels.jp <- c("9������", "�h�q�͋���", "���ۖ@���p�~", "�Ζk���N�Z�a", "�푈�Ӎߔ���", "��d���Ћ֎~", 
                     "��������", "�����ێ�", "���Y�p�~", "�����i�o����", "�v�w�ʐ�", "�Ƒ�����", "������", 
                     "�����p�~", "�i������")
png("Figure_1.png", 
    width = 4.2, height = 3.6, units = "in", pointsize = 7, res = 1000)
layout(matrix(1:15, 3, 5, byrow = TRUE))
par(mar = c(2, 2, 2, 1), lwd = 0.5)
for (i in 1:15) {
  plot(NULL, NULL, bty = "n", xlim = c(0.5, 5.5), ylim = c(0, 0.6), 
       main = issue.labels.jp[i], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  for (j in 1:5) {
    polygon(c(j - 0.5, j + 0.5, j + 0.5, j - 0.5), 
            c(0, 0, adult.issue.distribution[i, j], adult.issue.distribution[i, j]), 
            border = "#80808080", col = "#80808080")
    polygon(c(j - 0.5, j + 0.5, j + 0.5, j - 0.5), 
            c(0, 0, young.issue.distribution[i, j], young.issue.distribution[i, j]), lwd = 0.5)
  }
  text(5, 0.6, sprintf("%3.2f", weighted.mean(is.na(young.issue.data[, i]), weight = young.data$weight)))
  text(5, 0.54, sprintf("%3.2f", weighted.mean(is.na(adult.issue.data[, i]), weight = adult.data$weight)))
  axis(side = 1, at = c(1:5), cex = 0.8, lwd = 0.5)
  axis(side = 2, at = seq(0, 0.8, 0.2), cex = 0.8, lwd = 0.5)
}
dev.off()

## �{���ň���Ȃ��������_����
# ��N�w�����̖f�Վ��R���ƑS�N�����TPP�̗�ɂ̓_�~�[�f�[�^��}�����Ă���
young.issue.econ.data <- cbind(young.data[, 16:18], sample(1:5, nrow(young.data), replace = TRUE))
young.issue.econ.data[young.issue.econ.data == 6] <- NA
adult.issue.econ.data <- cbind(adult.data[, 16:17], sample(1:5, nrow(adult.data), replace = TRUE), adult.data[, 18])
adult.issue.econ.data[adult.issue.econ.data == 6] <- NA

issue.econ.distribution.young <- t(apply(young.issue.econ.data, 2, wpct, weight = young.data$weight))
issue.econ.distribution.adult <- t(apply(adult.issue.econ.data, 2, wpct, weight = adult.data$weight))

## �}A1
issue.econ.labels.jp <- c("�����", "��������", "TPP", "�f�Վ��R��")
cairo_pdf("Figure_A1.pdf", width = 5.04, height = 1.8, pointsize = 9, family = "Meiryo")
layout(matrix(1:4, 1, 4))
par(mar = c(2, 2, 2, 1), lwd = 0.5)
for (i in 1:4) {
  plot(NULL, NULL, bty = "n", xlim = c(0.5, 5.5), ylim = c(0, 0.6), 
       main = issue.econ.labels.jp[i], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  for (j in 1:5) {
    if (i != 3) {
      polygon(c(j - 0.5, j + 0.5, j + 0.5, j - 0.5), 
              c(0, 0, issue.econ.distribution.adult[i, j], issue.econ.distribution.adult[i, j]), 
              border = "#80808080", col = "#80808080")
    }
    if (i != 4) {
      polygon(c(j - 0.5, j + 0.5, j + 0.5, j - 0.5), 
              c(0, 0, issue.econ.distribution.young[i, j], issue.econ.distribution.young[i, j]), lwd = 0.5)
    }
  }
  if (i != 4) {
    text(5, 0.6, sprintf("%3.2f", weighted.mean(is.na(young.issue.econ.data[, i]), weight = young.data$weight)))
  }
  if (i != 3) {
    text(5, 0.54, sprintf("%3.2f", weighted.mean(is.na(adult.issue.econ.data[, i]), weight = adult.data$weight)))
  }
  axis(side = 1, at = c(1:5), cex = 0.8, lwd = 0.5)
  axis(side = 2, at = seq(0, 0.8, 0.2), cex = 0.8, lwd = 0.5)
}
dev.off()

#### ���_�ԓx�̃C�f�I���M�[�I��ѐ��̕��� ####
## ��ѐ��X�R�A�Ȃǂ̍쐬
# �e���_�Ɋւ���ԓx���E���Ȃ�+1�C�����Ȃ�-1�C�ӌ��Ȃ��Ȃ�0�����Z���Ă���
issue.agree.right <- c(1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0)  # ���]����
young.issue.sorted <- matrix(NA, nrow(young.issue.data.raw), 15)
adult.issue.sorted <- matrix(NA, nrow(adult.issue.data.raw), 15)
for (i in 1:15) {
  if (issue.agree.right[i] == 1) {
    young.issue.sorted[, i] <- ifelse(young.issue.data.raw[, i] < 3, 1, 
                                      ifelse(young.issue.data.raw[, i] == 3 | young.issue.data.raw[, i] == 6, 0, -1))
    adult.issue.sorted[, i] <- ifelse(adult.issue.data.raw[, i] < 3, 1, 
                                      ifelse(adult.issue.data.raw[, i] == 3 | adult.issue.data.raw[, i] == 6, 0, -1))
  } else {
    young.issue.sorted[, i] <- ifelse(young.issue.data.raw[, i] > 3, 1, 
                                      ifelse(young.issue.data.raw[, i] == 3 | young.issue.data.raw[, i] == 6, 0, -1))
    adult.issue.sorted[, i] <- ifelse(adult.issue.data.raw[, i] > 3, 1, 
                                      ifelse(adult.issue.data.raw[, i] == 3 | adult.issue.data.raw[, i] == 6, 0, -1))
  }
}

# ���Z���ꂽ���_�̐�Βl�i�g���f�[�^�ɂ���O�̈�ѐ��X�R�A�j
young.data$consistency <- abs(rowSums(young.issue.sorted))
adult.data$consistency <- abs(rowSums(adult.issue.sorted))

# �ӌ��Ȃ��ł��鑈�_�̐�
young.data$no.opinion <- rowSums(young.issue.sorted == 0)
adult.data$no.opinion <- rowSums(adult.issue.sorted == 0)

# ��т��Ă��Ȃ����_�̐�
young.data$inconsistency <- 15 - young.data$consistency - young.data$no.opinion
adult.data$inconsistency <- 15 - adult.data$consistency - adult.data$no.opinion

# �g���f�[�^���i�O�p�v���b�g�̍�}�̓s����C���сE�ӌ��Ȃ��E��т̏��Ƃ��Ă���j
young.data$y <- DR_data(cbind(young.data$inconsistency, 
                              young.data$no.opinion, 
                              young.data$consistency))
adult.data$y <- DR_data(cbind(adult.data$inconsistency, 
                              adult.data$no.opinion, 
                              adult.data$consistency))

## �f�B���N����A����
# ����
young.DR.result <- DirichReg(y ~ age + gender | age + gender, 
                             data = young.data, model = "alternative", weights = weight)
adult.DR.result <- DirichReg(y ~ age + gender | age + gender, 
                             data = adult.data, model = "alternative", weights = weight)

# �W���\�i�\A2�Cv2���ӌ��Ȃ���v3����тł��邱�Ƃɒ��ӂ��ꂽ���j
young.DR.result.summary <- summary(young.DR.result)
adult.DR.result.summary <- summary(adult.DR.result)

round(young.DR.result.summary$coef.mat, 3)
round(adult.DR.result.summary$coef.mat, 3)

## ���K���z�ߎ��ɂ��V�~�����[�V����
# �V�~�����[�V�����̂��߂̊֐�
DR.simulation <- function(result, newdata, iter = 1000, alpha = 0.05) {
  coef <- result$coefficients
  vcv <- result$vcov
  par.draw <- mvrnorm(iter, coef, vcv)
  coef.v2 <- par.draw[, stri_detect_regex(colnames(par.draw), "v2")]
  coef.v3 <- par.draw[, stri_detect_regex(colnames(par.draw), "v3")]
  exp.lp.v2 <- exp(coef.v2 %*% newdata)
  exp.lp.v3 <- exp(coef.v3 %*% newdata)
  mu.v1 <- 1 / (1 + exp.lp.v2 + exp.lp.v3)
  mu.v2 <- exp.lp.v2 / (1 + exp.lp.v2 + exp.lp.v3)
  mu.v3 <- exp.lp.v3 / (1 + exp.lp.v2 + exp.lp.v3)
  out <- cbind(mu.v1, mu.v2, mu.v3)
}

# 15�̑��_�ԓx���ڂ��Ɨ����Ɖ��肵���Ƃ��̈�ѐ��X�R�A�Ȃǂ̒l�̃V�~�����[�V����
young.issue.raw.distribution <- t(apply(young.issue.data.raw, 2, wpct, weight = young.data$weight))
sim.composition.young <- matrix(NA, 10000, 3)
for (i in 1:10000) {
  sim.responses <- rep(NA, 15)
  issue.sorted.sim <- rep(NA, 15)
  for (j in 1:15) {
    sim.response <- crossprod(1:6, rmultinom(1, 1, young.issue.raw.distribution[j, ]))
    if (issue.agree.right[j] == 1) {
      issue.sorted.sim[j] <- ifelse(sim.response < 3, 1, 
                                    ifelse(sim.response == 3 | sim.response == 6, 0, -1))
    } else {
      issue.sorted.sim[j] <- ifelse(sim.response > 3, 1, 
                                    ifelse(sim.response == 3 | sim.response == 6, 0, -1))
    }
  }
  sim.composition.young[i, 3] <- abs(sum(issue.sorted.sim)) / 15
  sim.composition.young[i, 2] <- sum(issue.sorted.sim == 0) / 15
  sim.composition.young[i, 1] <- 1 - sim.composition.young[i, 2] - sim.composition.young[i, 3]
}

adult.issue.raw.distribution <- t(apply(adult.issue.data.raw, 2, wpct, weight = adult.data$weight))
sim.composition.adult <- matrix(NA, 10000, 3)
for (i in 1:10000) {
  sim.responses <- rep(NA, 15)
  issue.sorted.sim <- rep(NA, 15)
  for (j in 1:15) {
    sim.response <- crossprod(1:6, rmultinom(1, 1, adult.issue.raw.distribution[j, ]))
    if (issue.agree.right[j] == 1) {
      issue.sorted.sim[j] <- ifelse(sim.response < 3, 1, 
                                    ifelse(sim.response == 3 | sim.response == 6, 0, -1))
    } else {
      issue.sorted.sim[j] <- ifelse(sim.response > 3, 1, 
                                    ifelse(sim.response == 3 | sim.response == 6, 0, -1))
    }
  }
  sim.composition.adult[i, 3] <- abs(sum(issue.sorted.sim)) / 15
  sim.composition.adult[i, 2] <- sum(issue.sorted.sim == 0) / 15
  sim.composition.adult[i, 1] <- 1 - sim.composition.adult[i, 2] - sim.composition.adult[i, 3]
}

# �}2
png("Figure_2.png", width = 6, height = 3, units = "in", pointsize = 6, res = 1000)
layout(matrix(1:2, 1, 2))
par(mar = c(1, 1, 1, 1), lwd = 0.5)
plot(rcomp(c(1 / 3, 1 / 3, 1 / 3)), pch = NA, labels = c("����", "�ӌ��Ȃ�", "���"), axes = TRUE)
isoPortionLines(labs = FALSE, lty = 3, col = "gray80")
for (i in 1:8) {
  sim.result <- rcomp(DR.simulation(young.DR.result, c(1, i + 14, 0)))
  ellipses(mean(sim.result), var(sim.result), sqrt(qchisq(0.95, 2)), col = "gray")
}
for (i in 1:8) {
  sim.result <- rcomp(DR.simulation(young.DR.result, c(1, i + 14, 0)))
  mn <- mean(sim.result)
  plot(mn, pch = 19, add = TRUE, labels = rep("", 3))
  if (i == 1) {
    text((2 * mn[2] + mn[3]) / 2, sqrt(3) * mn[3] / 2, "15M", pos = 2)
  }
  if (i == 8) {
    text((2 * mn[2] + mn[3]) / 2, sqrt(3) * mn[3] / 2, "22M", pos = 3)
  }
}
for (i in 1:8) {
  sim.result <- rcomp(DR.simulation(young.DR.result, c(1, i + 14, 1)))
  ellipses(mean(sim.result), var(sim.result), sqrt(qchisq(0.95, 2)), col = "gray")
}
for (i in 1:8) {
  sim.result <- rcomp(DR.simulation(young.DR.result, c(1, i + 14, 1)))
  mn <- mean(sim.result)
  plot(mn, pch = 19, add = TRUE, labels = rep("", 3))
  if (i == 1) {
    text((2 * mn[2] + mn[3]) / 2, sqrt(3) * mn[3] / 2, "15F", pos = 1)
  }
  if (i == 8) {
    text((2 * mn[2] + mn[3]) / 2, sqrt(3) * mn[3] / 2, "22F", pos = 4)
  }
}
plot(rcomp(colMeans(sim.composition.young)), pch = 4, add = TRUE, labels = rep("", 3))
text(0.5, 1, "��N�w����", font = 2, cex = 1.5)
plot(rcomp(c(1 / 3, 1 / 3, 1 / 3)), pch = NA, labels = c("����", "�ӌ��Ȃ�", "���"), axes = TRUE)
isoPortionLines(labs = FALSE, lty = 3, col = "gray80")
for (i in 1:7) {
  sim.result <- rcomp(DR.simulation(adult.DR.result, c(1, 10 * i + 9, 0)))
  ellipses(mean(sim.result), var(sim.result), sqrt(qchisq(0.95, 2)), col = "gray")
}
for (i in 1:7) {
  sim.result <- rcomp(DR.simulation(adult.DR.result, c(1, 10 * i + 9, 0)))
  mn <- mean(sim.result)
  plot(mn, pch = 19, add = TRUE, labels = rep("", 3))
  if (i == 1) {
    text((2 * mn[2] + mn[3]) / 2, sqrt(3) * mn[3] / 2, "19M", pos = 1)
  }
  if (i == 7) {
    text((2 * mn[2] + mn[3]) / 2, sqrt(3) * mn[3] / 2, "79M", pos = 3)
  }
}
for (i in 1:7) {
  sim.result <- rcomp(DR.simulation(adult.DR.result, c(1, 10 * i + 9, 1)))
  ellipses(mean(sim.result), var(sim.result), sqrt(qchisq(0.95, 2)), col = "gray")
}
for (i in 1:7) {
  sim.result <- rcomp(DR.simulation(adult.DR.result, c(1, 10 * i + 9, 1)))
  mn <- mean(sim.result)
  plot(mn, pch = 19, add = TRUE, labels = rep("", 3))
  if (i == 1) {
    text((2 * mn[2] + mn[3]) / 2, sqrt(3) * mn[3] / 2, "19F", pos = 1)
  }
  if (i == 7) {
    text((2 * mn[2] + mn[3]) / 2, sqrt(3) * mn[3] / 2, "79F", pos = 3)
  }
}
plot(rcomp(colMeans(sim.composition.adult)), pch = 4, add = TRUE, labels = rep("", 3))
text(0.5, 1, "�S�N���", font = 2, cex = 1.5)
dev.off()

## �N���2�����ꂽ����
young.age2.DR.result <- DirichReg(y ~ age + I(age ^ 2) + gender | age + I(age ^ 2) + gender, 
                                  data = young.data, model = "alternative", weights = weight)
adult.age2.DR.result <- DirichReg(y ~ age + I(age ^ 2) + gender | age + I(age ^ 2) + gender, 
                                  data = adult.data, model = "alternative", weights = weight)

young.age2.DR.result.summary <- summary(young.age2.DR.result)
adult.age2.DR.result.summary <- summary(adult.age2.DR.result)

# �\A3
round(young.age2.DR.result.summary$coef.mat, 3)
round(adult.age2.DR.result.summary$coef.mat, 3)

## �N��Ɛ��ʂ̌��ݍ�p������ꂽ����
young.int.DR.result <- DirichReg(y ~ age * gender | age * gender, 
                                 data = young.data, model = "alternative", weights = weight)
adult.int.DR.result <- DirichReg(y ~ age * gender | age * gender, 
                                 data = adult.data, model = "alternative", weights = weight)

young.int.DR.result.summary <- summary(young.int.DR.result)
adult.int.DR.result.summary <- summary(adult.int.DR.result)

# �\A4
round(young.int.DR.result.summary$coef.mat, 3)
round(adult.int.DR.result.summary$coef.mat, 3)

#### �C�f�I���M�[�I�I�D�ƒ����I�}�h���̈�v�x ####
## ���}�i�c���j�̗��z�_�̐���
# �����ƒ����f�[�^�̑��_�ԓx���ڂ̒��o
elite.response <- cbind(elite.data[, c("Q4_1", "Q4_3", "Q4_11", "Q4_14", "Q4_15", "Q5_3", "Q5_1")])
colnames(elite.response) <- c("defense", "North.Korea", "security", "surnames", "gay", "nuclear", "equality")

# ���]���ڂ̏���
elite.response$defense <- 6 - elite.response$defense
elite.response$North.Korea <- 6 - elite.response$North.Korea
elite.response$security <- 6 - elite.response$security
elite.response$equality <- 6 - elite.response$equality

# ���͑ΏۂƂȂ鍀�ڑS�Ă������l�̋c�������O����
elite.data.na.omit <- subset(elite.data, rowSums(is.na(elite.response)) < 7)
elite.response.na.omit <- subset(elite.response, rowSums(is.na(elite.response)) < 7)
nrow(elite.response.na.omit)

# �����t�����ڔ������_���f���̐���Ƃ��̌��ʂɊ�Â����z�_�̐���
IRT.result <- grm(elite.response.na.omit, IRT.param = FALSE)
elite.data.na.omit$theta <- factor.scores(IRT.result, resp.pattern = elite.response.na.omit)$score.dat$z1

# �����t�����ڔ������_���f���̐��茋�ʁi�\A6�j
round(do.call("rbind", IRT.result$coefficients), 3)

# �e�}�c���̗��z�_�̕��ϒl�i��22�j
LDP.theta <- mean(elite.data.na.omit$theta[elite.data.na.omit$PARTY == 1])  # �����}
DP.theta <- mean(elite.data.na.omit$theta[! is.na(elite.data.na.omit$MINSHIN)])  # ���i�}�i�����O�j
POH.theta <- mean(elite.data.na.omit$theta[elite.data.na.omit$PARTY == 7])  # ��]�̓}
CDP.theta <- mean(elite.data.na.omit$theta[elite.data.na.omit$PARTY == 8])  # ��������}
CGP.theta <- mean(elite.data.na.omit$theta[elite.data.na.omit$PARTY == 2])  # �����}
JCP.theta <- mean(elite.data.na.omit$theta[elite.data.na.omit$PARTY == 3])  # ���Y�}
JIP.theta <- mean(elite.data.na.omit$theta[elite.data.na.omit$PARTY == 4])  # ���{�ېV�̉�
LP.theta <- mean(elite.data.na.omit$theta[elite.data.na.omit$NAME %in% 
                                            c("�������ꐳ", "���򁁈�Y", "����������", 
                                              "�ʏ遁�f�j�[", "���ぁ�j�D", "���g���Y��")])  # ���R�}
SDP.theta <- mean(elite.data.na.omit$theta[elite.data.na.omit$PARTY == 5])  # �Ж��}

round(sort(data.frame(LDP.theta, DP.theta, POH.theta, CDP.theta, 
                      CGP.theta, JCP.theta, JIP.theta, LP.theta, SDP.theta)), 2)

## �I�������v���r�b�g���f��
# ��N�w�����ƑS�N����ŗ��p���鑈�_�ԓx���ڂ̒��o
young.response <- young.issue.data[, c("defense", "North.Korea", "security", "surnames", 
                                       "gay", "nuclear", "equality")]
adult.response <- adult.issue.data[, c("defense", "North.Korea", "security", "surnames", 
                                       "gay", "nuclear", "equality")]

# ���]���ڂ̏���
young.response$defense <- 6 - young.response$defense
young.response$security <- 6 - young.response$security
adult.response$defense <- 6 - adult.response$defense
adult.response$security <- 6 - adult.response$security

# ��N�w�����ƑS�N����ł̗��z�_����
young.data$theta <- factor.scores(IRT.result, resp.pattern = young.response)$score.dat$z1
adult.data$theta <- factor.scores(IRT.result, resp.pattern = adult.response)$score.dat$z1

# ���Ԗڂɋ߂����}�ɑ΂��ē}�h���������Ă��邩��\���ϐ��̍쐬
# �i���̎��_�ł́C���h�ɑ΂��ē}�h�������҂�}�h�����Ȃ��҂�0�ɂȂ�j
young.data$proximity.rank.raw <- rowSums(cbind(t(apply(cbind(abs(LDP.theta - young.data$theta), 
                                                             abs(DP.theta - young.data$theta), 
                                                             abs(CGP.theta - young.data$theta), 
                                                             abs(JCP.theta - young.data$theta), 
                                                             abs(JIP.theta - young.data$theta), 
                                                             abs(SDP.theta - young.data$theta), 
                                                             abs(LP.theta - young.data$theta)), 1, rank)), NA, NA, NA) * 
                                           diag(10)[young.data$partisanship, ], na.rm = TRUE)
adult.data$proximity.rank.raw <- rowSums(cbind(t(apply(cbind(abs(LDP.theta - adult.data$theta), 
                                                             abs(POH.theta - adult.data$theta), 
                                                             abs(CDP.theta - adult.data$theta), 
                                                             abs(CGP.theta - adult.data$theta), 
                                                             abs(JCP.theta - adult.data$theta), 
                                                             abs(JIP.theta - adult.data$theta), 
                                                             abs(SDP.theta - adult.data$theta)), 1, rank)), NA, NA, NA) * 
                                           diag(10)[adult.data$partisanship, ], na.rm = TRUE)

# �]���ϐ��̍쐬�i0�������l�Ƃ���j
young.data$proximity.rank <- young.data$proximity.rank.raw
young.data$proximity.rank[young.data$proximity.rank == 0] <- NA
young.data$proximity.rank <- factor(young.data$proximity.rank)
adult.data$proximity.rank <- adult.data$proximity.rank.raw
adult.data$proximity.rank[adult.data$proximity.rank == 0] <- NA
adult.data$proximity.rank <- factor(adult.data$proximity.rank)

# �}�h���������Ƃ�\���ϐ��̍쐬
young.data$not.indep <- (young.data$proximity.rank.raw != 0) * 1
adult.data$not.indep <- (adult.data$proximity.rank.raw != 0) * 1

# ���_�ԓx���ڂɂ�����u�킩��Ȃ��v�̐��̍쐬
young.data$DK <- rowSums(young.issue.data.raw == 6)
adult.data$DK <- rowSums(adult.issue.data.raw == 6)

# �����l�ݒ�̂��߂̕��͂�survey.design�I�u�W�F�N�g�̍쐬
young.survey.design <- svydesign(id = ~ 1, weights = young.data$weight, data = young.data)
adult.survey.design <- svydesign(id = ~ 1, weights = adult.data$weight, data = adult.data)

# �ޓx�֐�
selection.ord.probit <- function(est,  # �p�����[�^
                                 y,  # �]���ϐ�
                                 x1,  # y.star_1�̓Ɨ��ϐ�
                                 x2,  # y.star_2�̓Ɨ��ϐ�
                                 w  # �E�F�C�g
) {
  n <- length(y)
  k <- max(y, na.rm = TRUE)
  y.dummy <- diag(k + 1)[y + 1, ]
  beta1 <- est[1:ncol(x1)]
  beta2 <- est[(ncol(x1) + 1):(ncol(x1) + ncol(x2))]
  tau <- c(0, rep(NA, k - 2))
  for (i in 2:(k - 1)) {
    tau[i] <- tau[i - 1] + exp(est[ncol(x1) + ncol(x2) + i - 1])
  }
  rho <- 2 * pnorm(est[length(est)]) - 1
  beta1.x1 <- c(x1 %*% beta1)
  beta2.x2 <- c(x2 %*% beta2)
  probs <- matrix(NA, n, k + 1)
  probs[, 1] <- pnorm(-1 * beta1.x1)
  probs[, 2] <- pbivnorm(beta1.x1, -1 * beta2.x2, rho = -1 * rho)
  for (i in 1:(k - 2)) {
    probs[, i + 2] <- pbivnorm(beta1.x1, tau[i + 1] - beta2.x2, rho = -1 * rho) - 
      pbivnorm(beta1.x1, tau[i] - beta2.x2, rho = -1 * rho)
  }
  probs[, k + 1] <- pbivnorm(beta1.x1, beta2.x2 - tau[k - 1], rho = rho)
  ll <- w * rowSums(y.dummy * log(probs))
  return(ll)
}

# �Ŗސ���̂��߂̊֐�
SOP.analysis <- function(y.var,  # �]���ϐ��̕ϐ���
                         formula.1,  # y.star_1�̎�
                         formula.2,  # y.star_2�̎�
                         design,  # ���͂���survey.design�I�u�W�F�N�g
                         iter  # ���s��
) {
  y <- design$variables[, y.var]
  y.numeric <- as.numeric(y)
  x1 <- model.matrix(formula.1, design$variables)
  x2 <- model.matrix(formula(paste(as.character(formula.1)[2], "~", as.character(formula.2)[3])), 
                     design$variables)
  # y.star_1�̎��̃p�����[�^�̏����l�ݒ�
  bp.result <- svyglm(formula.1, design = design, 
                      family = quasibinomial(link = "probit"))
  # y.star_2�̎��̃p�����[�^�̏����l�ݒ�
  op.result <- svyolr(formula.2, design = design, 
                      method = "probit")
  base.inits <- c(bp.result$coefficients, 
                  -1 * op.result$zeta[1], op.result$coefficients, 
                  rep(NA, 5), 0)
  for (i in 1:5) {  # 臒l�p�����[�^�̎��ʕ��@��ς��邽�߂̕ϊ�
    base.inits[i + ncol(x1) + ncol(x2)] <- log(op.result$zeta[i + 1] - op.result$zeta[i])
  }
  # �Ŗސ���
  MLE.result.list <- list()
  MLE.par.matrix <- matrix(NA, iter, length(base.inits) + 1)
  for (i in 1:iter) {
    set.seed(1000 * i)
    if (iter == 1) {
      inits <- base.inits
    } else {  # iter > 1�̏ꍇ�͏����l��ς��ČJ��Ԃ�
      inits <- base.inits + rnorm(ncol(x1) + ncol(x2) + 6, 0, 0.1)
    }
    MLE.result.list[[i]] <- maxBFGS(fn = selection.ord.probit, start = inits, 
                                    y = y.numeric, x1 = x1, x2 = x2, w = design$prob ^ -1)
    if (MLE.result.list[[i]]$message == "successful convergence ") {
      MLE.par.matrix[i, 1:(ncol(x1) + ncol(x2))] <- MLE.result.list[[i]]$estimate[1:(ncol(x1) + ncol(x2))]
      tau <- c(0, rep(NA, 5))
      for (j in 2:6) {
        tau[j] <- tau[j - 1] + exp(MLE.result.list[[i]]$estimate[j + ncol(x1) + ncol(x2) - 1])
      }
      MLE.par.matrix[i, (ncol(x1) + ncol(x2) + 1):(ncol(x1) + ncol(x2) + 5)] <- tau[2:6]
      MLE.par.matrix[i, ncol(x1) + ncol(x2) + 6] <- 2 * pnorm(MLE.result.list[[i]]$estimate[ncol(x1) + ncol(x2) + 6]) - 1
      MLE.par.matrix[i, ncol(x1) + ncol(x2) + 7] <- MLE.result.list[[i]]$maximum
    }
    cat("Iteration", i, "finished at", date(), "\n")
  }
  # iter > 1�̏ꍇ�͍ł��ޓx���傫�����̂��̗p����
  max.loglik <- which.max(MLE.par.matrix[, ncol(MLE.par.matrix)])
  # �_����l
  par <- MLE.result.list[[max.loglik]]$estimate
  # �d�ݕt�����U�����U�s��
  vcov <- solve(MLE.result.list[[max.loglik]]$hessian) %*% 
    (crossprod(MLE.result.list[[max.loglik]]$gradientObs)) %*% 
    solve(MLE.result.list[[max.loglik]]$hessian)
  # �W���덷
  se <- sqrt(diag(vcov))
  # 臒l�p�����[�^�̃V�~�����[�V����
  set.seed(12345)
  par.draw <- mvrnorm(1000, par, vcov)
  tau <- matrix(NA, 1000, 6)
  tau[, 1] <- 0
  for (i in 2:6) {
    tau[, i] <- tau[, i - 1] + exp(par.draw[, ncol(x1) + ncol(x2) + i - 1])
  }
  # ���փp�����[�^�̃V�~�����[�V����
  rho <- 2 * pnorm(par.draw[, ncol(x1) + ncol(x2) + 6]) - 1
  result <- list(summary = list(coef.1 = par[1:ncol(x1)], coef.1.se = se[1:ncol(x1)], 
                                coef.1.p = 2 * (1 - pnorm(abs(par[1:ncol(x1)]) / se[1:ncol(x1)])), 
                                coef.2 = par[(ncol(x1) + 1):(ncol(x1) + ncol(x2))], 
                                coef.2.se = se[(ncol(x1) + 1):(ncol(x1) + ncol(x2))], 
                                coef.2.p = 2 * (1 - pnorm(abs(par[(ncol(x1) + 1):(ncol(x1) + ncol(x2))]) / 
                                                            se[(ncol(x1) + 1):(ncol(x1) + ncol(x2))])), 
                                tau = colMeans(tau), tau.se = apply(tau, 2, sd), 
                                rho = mean(rho), rho.se = sd(rho)), 
                 par.matrix = MLE.par.matrix, par = par, vcov = vcov)
}

# ����
young.SOP.result <- SOP.analysis(y.var = "proximity.rank.raw", 
                                 formula.1 = not.indep ~ age + gender + DK, 
                                 formula.2 = proximity.rank ~ age + gender, 
                                 design = young.survey.design, iter = 10)

adult.SOP.result <- SOP.analysis(y.var = "proximity.rank.raw", 
                                 formula.1 = not.indep ~ age + gender + DK, 
                                 formula.2 = proximity.rank ~ age + gender, 
                                 design = adult.survey.design, iter = 10)

# �W���\�i�\A7�j
lapply(young.SOP.result$summary, function(x) round(x, 3))
lapply(adult.SOP.result$summary, function(x) round(x, 3))

# �����l�ɂ���Č��ʂ��傫���قȂ�Ȃ����Ƃ̊m�F
round(young.SOP.result$par.matrix, 3)
round(adult.SOP.result$par.matrix, 3)

## ���K���z�ߎ��ɂ��V�~�����[�V����
# �V�~�����[�V�����̂��߂̊֐�
SOP.simulation <- function(par,  # �p�����[�^
                           vcv,  # ���U�����U�s��
                           new.x1,  # y.star_1�̓Ɨ��ϐ�
                           new.x2,  # y.star_2�̓Ɨ��ϐ�
                           pred,  # Pr(y = pred)�܂���Pr(y < pred)�����߂�
                           iter = 1000,  # ���s��
                           alpha = 0.05,  # �M����Ԃ̐���
                           seed = 12345,  # �����̎�
                           conditional = FALSE,  # y.star_1 >= 0�������t���邩�ۂ�
                           cumulative = FALSE  # Pr(y < pred)�Ȃ�TRUE�CPr(y = pred)�Ȃ�FALSE
) {
  set.seed(seed)
  # �V�����p�����[�^���T���v�����O����
  par.draw <- mvrnorm(iter, par, vcv)
  beta1 <- par.draw[, 1:length(new.x1)]
  beta2 <- par.draw[, (length(new.x1) + 1):(length(new.x1) + length(new.x2))]
  tau <- matrix(NA, iter, 6)
  tau[, 1] <- 0
  for (i in 2:6) {
    tau[, i] <- tau[, i - 1] + exp(par.draw[, length(new.x1) + length(new.x2) + i - 1])
  }
  rho <- 2 * pnorm(par.draw[, length(new.x1) + length(new.x2) + 6]) - 1
  # �\���m�������߂�
  beta1.x1 <- c(beta1 %*% new.x1)
  beta2.x2 <- c(beta2 %*% new.x2)
  probs <- matrix(NA, iter, 8)
  probs[, 1] <- pnorm(-1 * beta1.x1)
  probs[, 2] <- pbivnorm(beta1.x1, -1 * beta2.x2, rho = -1 * rho)
  for (j in 1:5) {
    probs[, j + 2] <- pbivnorm(beta1.x1, tau[, j + 1] - beta2.x2, rho = -1 * rho) - 
      pbivnorm(beta1.x1, tau[, j] - beta2.x2, rho = -1 * rho)
  }
  probs[, 8] <- pbivnorm(beta1.x1, beta2.x2 - tau[, 6], rho = rho)
  if (conditional == TRUE) {  # y.star_1 >= 0�ŏ����t�����m���̏o��
    if (cumulative == TRUE) {
      conditional.probs <- rowSums(probs[, 2:(pred + 1)]) / (1 - probs[, 1])
    } else {
      conditional.probs <- probs[, pred + 1] / (1 - probs[, 1])
    }
    out <- matrix(NA, 1, 3)
    out[, 1] <- mean(conditional.probs)
    out[, 2:3] <- quantile(conditional.probs, probs = c(alpha / 2, 1 - alpha / 2))
  } else {  # �����t���łȂ��m���̏o��
    out <- matrix(NA, 1, 3)
    if (cumulative == TRUE) {
      out[, 1] <- mean(rowSums(probs[, 2:(pred + 1)]))
      out[, 2:3] <- quantile(rowSums(probs[, 2:(pred + 1)]), probs = c(alpha / 2, 1 - alpha / 2))
    } else {
      out[, 1] <- mean(probs[, pred + 1])
      out[, 2:3] <- quantile(probs[, pred + 1], probs = c(alpha / 2, 1 - alpha / 2))
    }
  }
  return(out)
}

# �V�~�����[�V�����̎��s
SOP.sim.plot.young <- array(NA, c(8, 3, 3, 2))
for (i in 1:8) {  # �N���1�΂��ω�������
  # ��v7���}�̂����ꂩ�ɓ}�h�������m���E�j��
  SOP.sim.plot.young[i, 1, , 1] <- SOP.simulation(young.SOP.result$par, young.SOP.result$vcov, 
                                                  c(1, i + 14, 0, 1), c(1, i + 14, 0), 7, 
                                                  cumulative = TRUE)
  # ��v7���}�̂����ꂩ�ɓ}�h�������m���E����
  SOP.sim.plot.young[i, 1, , 2] <- SOP.simulation(young.SOP.result$par, young.SOP.result$vcov, 
                                                  c(1, i + 14, 1, 1), c(1, i + 14, 1), 7, 
                                                  cumulative = TRUE)
  # �ł��߂����}�ɑ΂��ē}�h�����������t���m���E�j��
  SOP.sim.plot.young[i, 2, , 1] <- SOP.simulation(young.SOP.result$par, young.SOP.result$vcov, 
                                                  c(1, i + 14, 0, 1), c(1, i + 14, 0), 1, 
                                                  conditional = TRUE)
  # �ł��߂����}�ɑ΂��ē}�h�����������t���m���E����
  SOP.sim.plot.young[i, 2, , 2] <- SOP.simulation(young.SOP.result$par, young.SOP.result$vcov, 
                                                  c(1, i + 14, 1, 1), c(1, i + 14, 1), 1, 
                                                  conditional = TRUE)
  # 3�Ԗڈȓ��ɋ߂����}�ɑ΂��ē}�h�����������t���m���E�j��
  SOP.sim.plot.young[i, 3, , 1] <- SOP.simulation(young.SOP.result$par, young.SOP.result$vcov, 
                                                  c(1, i + 14, 0, 1), c(1, i + 14, 0), 3, 
                                                  conditional = TRUE, cumulative = TRUE)
  # 3�Ԗڈȓ��ɋ߂����}�ɑ΂��ē}�h�����������t���m���E����
  SOP.sim.plot.young[i, 3, , 2] <- SOP.simulation(young.SOP.result$par, young.SOP.result$vcov, 
                                                  c(1, i + 14, 1, 1), c(1, i + 14, 1), 3, 
                                                  conditional = TRUE, cumulative = TRUE)
}

SOP.sim.plot.adult <- array(NA, c(62, 3, 3, 2))
for (i in 1:62) {  # �N���1�΂��ω�������
  # ��v7���}�̂����ꂩ�ɓ}�h�������m���E�j��
  SOP.sim.plot.adult[i, 1, , 1] <- SOP.simulation(adult.SOP.result$par, adult.SOP.result$vcov, 
                                                  c(1, i + 17, 0, 1), c(1, i + 17, 0), 7, 
                                                  cumulative = TRUE)
  # ��v7���}�̂����ꂩ�ɓ}�h�������m���E����
  SOP.sim.plot.adult[i, 1, , 2] <- SOP.simulation(adult.SOP.result$par, adult.SOP.result$vcov, 
                                                  c(1, i + 17, 1, 1), c(1, i + 17, 1), 7, 
                                                  cumulative = TRUE)
  # �ł��߂����}�ɑ΂��ē}�h�����������t���m���E�j��
  SOP.sim.plot.adult[i, 2, , 1] <- SOP.simulation(adult.SOP.result$par, adult.SOP.result$vcov, 
                                                  c(1, i + 17, 0, 1), c(1, i + 17, 0), 1, 
                                                  conditional = TRUE)
  # �ł��߂����}�ɑ΂��ē}�h�����������t���m���E����
  SOP.sim.plot.adult[i, 2, , 2] <- SOP.simulation(adult.SOP.result$par, adult.SOP.result$vcov, 
                                                  c(1, i + 17, 1, 1), c(1, i + 17, 1), 1, 
                                                  conditional = TRUE)
  # 3�Ԗڈȓ��ɋ߂����}�ɑ΂��ē}�h�����������t���m���E�j��
  SOP.sim.plot.adult[i, 3, , 1] <- SOP.simulation(adult.SOP.result$par, adult.SOP.result$vcov, 
                                                  c(1, i + 17, 0, 1), c(1, i + 17, 0), 3, 
                                                  conditional = TRUE, cumulative = TRUE)
  # 3�Ԗڈȓ��ɋ߂����}�ɑ΂��ē}�h�����������t���m���E����
  SOP.sim.plot.adult[i, 3, , 2] <- SOP.simulation(adult.SOP.result$par, adult.SOP.result$vcov, 
                                                  c(1, i + 17, 1, 1), c(1, i + 17, 1), 3, 
                                                  conditional = TRUE, cumulative = TRUE)
}

# �}3
png("Figure_3_lower.png", width = 5, height = 2, units = "in", pointsize = 7, res = 1000)
layout(matrix(1:3, 1, 3))
par(mar = c(4, 4, 2.5, 1), oma = c(0, 0, 3, 0), lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(15, 22), ylim = c(0, 1), 
     main = "��v7���}�̂����ꂩ�ɑ΂���\n�}�h��������", xlab = "�N��", ylab = "�m��", 
     xaxt = "n", yaxt = "n", cex.main = 1)
polygon(c(15:22, 22:15), 
        c(SOP.sim.plot.young[, 1, 2, 1], rev(SOP.sim.plot.young[, 1, 3, 1])), 
        border = NA, col = "#A0A0A080")
polygon(c(15:22, 22:15), 
        c(SOP.sim.plot.young[, 1, 2, 2], rev(SOP.sim.plot.young[, 1, 3, 2])), 
        border = NA, col = "#A0A0A080")
lines(15:22, SOP.sim.plot.young[, 1, 1, 1])
lines(15:22, SOP.sim.plot.young[, 1, 1, 2], lty = 2)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(15, 22), ylim = c(0, 1), 
     main = "�ł��߂����}�ɑ΂���\n�}�h��������", xlab = "�N��", ylab = "�m��", 
     xaxt = "n", yaxt = "n", cex.main = 1)
polygon(c(15:22, 22:15), 
        c(SOP.sim.plot.young[, 2, 2, 1], rev(SOP.sim.plot.young[, 2, 3, 1])), 
        border = NA, col = "#A0A0A080")
polygon(c(15:22, 22:15), 
        c(SOP.sim.plot.young[, 2, 2, 2], rev(SOP.sim.plot.young[, 2, 3, 2])), 
        border = NA, col = "#A0A0A080")
lines(15:22, SOP.sim.plot.young[, 2, 1, 1])
lines(15:22, SOP.sim.plot.young[, 2, 1, 2], lty = 2)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(15, 22), ylim = c(0, 1), 
     main = "3�Ԗڈȓ��ɋ߂����}�ɑ΂���\n�}�h��������", xlab = "�N��", ylab = "�m��", 
     xaxt = "n", yaxt = "n", cex.main = 1)
polygon(c(15:22, 22:15), 
        c(SOP.sim.plot.young[, 3, 2, 1], rev(SOP.sim.plot.young[, 3, 3, 1])), 
        border = NA, col = "#A0A0A080")
polygon(c(15:22, 22:15), 
        c(SOP.sim.plot.young[, 3, 2, 2], rev(SOP.sim.plot.young[, 3, 3, 2])), 
        border = NA, col = "#A0A0A080")
lines(15:22, SOP.sim.plot.young[, 3, 1, 1])
lines(15:22, SOP.sim.plot.young[, 3, 1, 2], lty = 2)
title("��N�w����", outer = TRUE, line = 1, cex = 1.2)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
dev.off()

png("Figure_3_upper.png", width = 5, height = 2, units = "in", pointsize = 7, res = 1000)
layout(matrix(1:3, 1, 3))
par(mar = c(4, 4, 2.5, 1), oma = c(0, 0, 3, 0), lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(18, 79), ylim = c(0, 1), 
     main = "��v7���}�̂����ꂩ�ɑ΂���\n�}�h��������", xlab = "�N��", ylab = "�m��", 
     xaxt = "n", yaxt = "n", cex.main = 1)
polygon(c(18:79, 79:18), 
        c(SOP.sim.plot.adult[, 1, 2, 1], rev(SOP.sim.plot.adult[, 1, 3, 1])), 
        border = NA, col = "#A0A0A080")
polygon(c(18:79, 79:18), 
        c(SOP.sim.plot.adult[, 1, 2, 2], rev(SOP.sim.plot.adult[, 1, 3, 2])), 
        border = NA, col = "#A0A0A080")
lines(18:79, SOP.sim.plot.adult[, 1, 1, 1])
lines(18:79, SOP.sim.plot.adult[, 1, 1, 2], lty = 2)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(18, 79), ylim = c(0, 1), 
     main = "�ł��߂����}�ɑ΂���\n�}�h��������", xlab = "�N��", ylab = "�m��", 
     xaxt = "n", yaxt = "n", cex.main = 1)
polygon(c(18:79, 79:18), 
        c(SOP.sim.plot.adult[, 2, 2, 1], rev(SOP.sim.plot.adult[, 2, 3, 1])), 
        border = NA, col = "#A0A0A080")
polygon(c(18:79, 79:18), 
        c(SOP.sim.plot.adult[, 2, 2, 2], rev(SOP.sim.plot.adult[, 2, 3, 2])), 
        border = NA, col = "#A0A0A080")
lines(18:79, SOP.sim.plot.adult[, 2, 1, 1])
lines(18:79, SOP.sim.plot.adult[, 2, 1, 2], lty = 2)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", xlim = c(18, 79), ylim = c(0, 1), 
     main = "3�Ԗڈȓ��ɋ߂����}�ɑ΂���\n�}�h��������", xlab = "�N��", ylab = "�m��", 
     xaxt = "n", yaxt = "n", cex.main = 1)
polygon(c(18:79, 79:18), 
        c(SOP.sim.plot.adult[, 3, 2, 1], rev(SOP.sim.plot.adult[, 3, 3, 1])), 
        border = NA, col = "#A0A0A080")
polygon(c(18:79, 79:18), 
        c(SOP.sim.plot.adult[, 3, 2, 2], rev(SOP.sim.plot.adult[, 3, 3, 2])), 
        border = NA, col = "#A0A0A080")
lines(18:79, SOP.sim.plot.adult[, 3, 1, 1])
lines(18:79, SOP.sim.plot.adult[, 3, 1, 2], lty = 2)
title("�S�N���", outer = TRUE, line = 1, cex = 1.2)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
dev.off()

## y.star_2��DK���������ւł��邱�Ƃ̊m�F
# �P���ȃv���r�b�g���f���ɂ�����DK���̌W���̊m�F�i�\A8�j
young.DK.OP.result <- svyolr(proximity.rank ~ age + gender + DK, design = young.survey.design, 
                             method = "probit")
young.DK.OP.result.summary <- summary(young.DK.OP.result)
round(young.DK.OP.result.summary$coefficients, 3)
round((1 - pnorm(abs(young.DK.OP.result.summary$coefficients[, 3]))) * 2, 3)

adult.DK.OP.result <- svyolr(proximity.rank ~ age + gender + DK, design = adult.survey.design, 
                             method = "probit")
adult.DK.OP.result.summary <- summary(adult.DK.OP.result)
round(adult.DK.OP.result.summary$coefficients, 3)
round((1 - pnorm(abs(adult.DK.OP.result.summary$coefficients[, 3]))) * 2, 3)

# �N���DK���̑��ցi��25�j
round(wtd.cor(young.data$age, young.data$DK, young.data$weight), 2)
round(wtd.cor(adult.data$age, adult.data$DK, adult.data$weight), 2)

## �u���\�E�v�̗p��ɑ΂���m���̎��ȔF������ꂽ���́i�\A9�j
young.knowledge.SOP.result <- SOP.analysis(y.var = "proximity.rank.raw", 
                                           formula.1 = not.indep ~ age + gender + knowledge + DK, 
                                           formula.2 = proximity.rank ~ age + gender + knowledge, 
                                           design = young.survey.design, iter = 1)
lapply(young.knowledge.SOP.result$summary, function(x) round(x, 3))

adult.knowledge.SOP.result <- SOP.analysis(y.var = "proximity.rank.raw", 
                                           formula.1 = not.indep ~ age + gender + knowledge + DK, 
                                           formula.2 = proximity.rank ~ age + gender + knowledge, 
                                           design = adult.survey.design, iter = 1)
lapply(adult.knowledge.SOP.result$summary, function(x) round(x, 3))

## �N���2�����ꂽ���́i�\A10�j
young.age2.SOP.result <- SOP.analysis(y.var = "proximity.rank.raw", 
                                      formula.1 = not.indep ~ age + I(age ^ 2) + gender + DK, 
                                      formula.2 = proximity.rank ~ age + I(age ^ 2) + gender, 
                                      design = young.survey.design, iter = 1)
lapply(young.age2.SOP.result$summary, function(x) round(x, 3))

adult.age2.SOP.result <- SOP.analysis(y.var = "proximity.rank.raw", 
                                      formula.1 = not.indep ~ age + I(age ^ 2) + gender + DK, 
                                      formula.2 = proximity.rank ~ age + I(age ^ 2) + gender, 
                                      design = adult.survey.design, iter = 1)
lapply(adult.age2.SOP.result$summary, function(x) round(x, 3))

## �N��Ɛ��ʂ̌��ݍ�p������ꂽ���́i�\A11�j
young.int.SOP.result <- SOP.analysis(y.var = "proximity.rank.raw", 
                                     formula.1 = not.indep ~ age * gender + DK, 
                                     formula.2 = proximity.rank ~ age * gender, 
                                     design = young.survey.design, iter = 1)
lapply(young.int.SOP.result$summary, function(x) round(x, 3))

adult.int.SOP.result <- SOP.analysis(y.var = "proximity.rank.raw", 
                                     formula.1 = not.indep ~ age * gender + DK, 
                                     formula.2 = proximity.rank ~ age * gender, 
                                     design = adult.survey.design, iter = 1)
lapply(adult.int.SOP.result$summary, function(x) round(x, 3))